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This  work  presents  a  three-dimensional,  steady-state,  non-isothermal  model  of  a  high-temperature 
polymer-electrolyte-membrane  fuel  cell  (HTPEMFC)  using  a  phosphoric  acid-doped  polybenzimidazole 
(PBI/H3PO4)  sol-gel  membrane.  The  model  accounts  for  the  gold-plated  copper  current  collector  plates, 
the  bipolar  plates,  all  gas  flow  channels  (flow-field),  the  gas  diffusion  layers,  the  reaction  layers,  and 
the  membrane.  Electrochemical  reactions  are  modeled  using  an  agglomerate  approach  and  include  the 
gas  diffusivity  and  the  gas  solubility.  The  conductivity  of  the  membrane  is  modeled  using  the  Arrhenius 
equation  to  describe  the  temperature  dependence.  Finite  elements  are  used  to  discretize  all  computa¬ 
tional  subdomains,  and  a  commercially  available  code  is  used  to  solve  the  problem.  The  predicted  values 
are  compared  to  typical  operating  conditions,  and  a  good  agreement  is  found.  The  current  density,  the 
solid-  and  fluid-(gas)-phase  temperatures  and  other  quantities  are  analyzed  throughout  the  computa¬ 
tional  subdomains.  It  was  observed  that  the  Arrhenius  approach  is  valid  in  a  certain  temperature  range 
and  may  overpredict  the  PBI/H3PO4  sol-gel  membrane  conductivity  at  higher  solid-phase  temperatures. 
Moreover,  it  is  shown  how  the  fluid-(gas)-phase  temperature  influences  the  solid-phase  temperature 
and  the  current  density  distribution.  Concrete  values  are  deduced  from  the  simulations  and  discussed 
according  to  experimental  test. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Within  the  last  four  decades,  tremendous  effort  has  been 
invested  in  both  theoretical  and  experimental  investigations 
of  low-temperature  polymer-electrolyte-membrane  fuel  cells 
(LTPEMFC).  Based  on  phosphoric  acid  fuel  cell  (PAFC)  technology, 
the  HTPEMFC  was  presented  as  an  interesting  alternative  in  appli¬ 
cations  such  as  small  stationary  power  plants  (combined  heat  and 
power,  CHP  or  piCHP)  and  power  trains.  Because  of  the  benefits  of 
using  higher  operating  temperatures  (e.g.,  160-1 80  °C),  this  tech¬ 
nology  may  play  a  significant  role  in  new  fuel  cell  systems  because 
carbon  monoxide  (CO)-rich  gases  in  the  range  of  several  percent¬ 
ages  can  be  directly  fed  into  the  cell. 

The  possibility  of  using  a  PBI/H3PO4  system  as  a  proton  conduc¬ 
tor  in  fuel  cells  was  introduced  by  Wainright  et  al.  [1  ],  Samms  et  al. 
[2],  and  Wang  et  al.  [3],  among  others.  They  tested  polybenzimi¬ 
dazole  films  doped  with  phosphoric  acid  (PBI/H3PO4)  as  potential 
polymer  electrolytes  for  fuel  cell  applications  under  various  oper¬ 
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ating  conditions.  Since  then,  the  physiochemical  properties,  water 
uptake,  conductivity  mechanisms  and  membrane  performance  has 
been  studied  by  several  groups  [4-7].  Xiao  et  al.  [8]  presented 
sol-gel  membranes  with  high  levels  of  phosphoric  acid,  high 
conductivities,  and  acceptable  mechanical  properties  at  elevated 
temperatures.  To  date,  only  a  few  publications  have  investigated 
the  short-  and  long-term  behavioral  patterns  of  commercially  avail¬ 
able  HTPEM  membrane-electrode-assemblies  (MEAs).  Schmidt  [9] 
summarized  HTPEMFC  durability  and  degradation  using  a  Celtec®- 
P  Series  1000  MEA.  The  properties  of  this  product  running  in 
start/stop  operation  mode  were  discussed  in  [10].  Electrochemi¬ 
cal  impedance  spectroscopy  (EIS)  was  used  by  Jalani  et  al.  [11]  to 
obtain  a  detailed  view  of  PBI/H3PO4  sol-gel  membrane  processes 
under  different  operating  conditions.  They  also  highlighted  the  use 
of  oxygen  as  cathode  gas  instead  of  air.  In  [12],  Zhang  et  al.  used  EIS, 
cyclic  voltammetry  (CV)  and  fuel  cell  performance  simulations  to 
obtain  the  exchange  current  densities  and  activation  energies  for 
both  half  cell  reactions.  Moreover,  they  proposed  different  methods 
to  improve  gas  diffusion  processes.  Zhang  et  al.  [  1 3  ]  and  Li  et  al.  [  1 4] 
summarized  the  current  status  of  HTPEMFC  research  and  develop¬ 
ment  and  described  the  advantages  and  challenges  of  operating  a 
fuel  cell  at  high  temperatures. 

Only  a  handful  of  publications  concerning  HTPEMFC  model¬ 
ing  are  currently  available.  Cheddie  et  al.  [15-18]  published  one-, 
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Nomenclature 

a  surface  area  (m2);  activity 

c  concentration  (mol  m-3) 

Cp  heat  capacity  (J  kg-1  K-1 ) 

d  average  hopping  distance  (m) 

D  diffusion  coefficient  (m2  s-1 ) 

E  potential  (V) 

A E  activation  energy  (J  mol-1 ) 

/  ratio  factor 

F  Faraday  constant  (A  s  mol-1 ) 

h  heat  transfer  coefficient  (W  m-3  K-1 )  (W  m-2  K-1 ) 

H  solubility  (mol  m-3  atm-1)  (mol  m-3  Pa-1) 

(molm-3  bar-1) 

i  current  density  (A  m-2 ) 

j  current  density  (volumetric)  (A  m-3 ) 

/  identity 

k  thermal  conductivity  (W m-1 1<-1 );  constant 

kp  permeability  (m2) 

I<  used  in  Eqs.  (7)  and  (8) 

l  thickness  (m) 

m  loading  (kgm-2);  mass  fraction 

M  molar  mass  (kg  mol-1 ) 

P  pressure  (Pa) 

r  radius  (m) 

R  gas  constant  (J  mol-1  K-1 ) 

S  general  source/sink  term  (kgm-3s-1)  (Am-3) 

(Wm-3) 

AS  entropy  (J  mol-1  K-1 ) 

T  phase  temperature  (°C)  (K) 

u  velocity  vector  (ms-1) 

x  mole  fraction;  dimension  in  x-direction  (m) 

X  membrane  doping  level 

y  dimension  in  y-direction  (m) 

z  dimension  in  z-direction  (m);  charge  number 

Greek  letters 

a  transfer  coefficient;  reciprocal  of  the  number  of  all 

possible  hopping  directions 
s  porosity 

0  phase  potential  (V) 

r\  dynamic  viscosity  (Pa  s);  overpotential  (V) 

v  factor  in  error  estimation;  diffusive  volume 

(m3  mol-1) 

v0  hopping  frequency 

p  density  (kgm-3) 

cr  conductivity  (Sm-1);  pre-exponential  factor 

(SKm-1)  (SKcm-1) 
co  mass  fraction 

Superscripts 

0  void,  reference 

a  activation 

eff  effective 

/  formation 

in  inlet 

OCV  open  circuit  voltage 

S  solid  fraction  in  GDL 

T  transposed 

Subscripts 

agg  agglomerate 

a  anode-side 

BPP  bipolar  plate 

c  cathode-side 


c 

carbon 

Cu 

copper 

/ 

fluid-phase 

GDL 

gas  diffusion  layer 

h2 

hydrogen 

h2o 

water 

H3PO4 

phosphoric  acid 

i 

index 

j 

index 

m 

membrane  phase 

02 

oxygen 

PBI 

polybenzimidazole 

PBI /(X- 

2  )H3  P04  PBI/amorphous  phase  FI3  P04 

Pt 

platinum 

Pt/C 

platinum  to  carbon  ratio 

RL 

reaction  layer 

s 

solid-phase 

Ts/Tf 

solid  to  fluid  heat  transfer 

two-,  and  three-dimensional  models  accounting  for  different  oper¬ 
ating  conditions,  membrane  and  reaction  layer  properties,  layout 
optimization,  and  gas  solubility  and  gaseous  dissolution  into 
the  aqueous  phase.  Steady-state  and  transient  three-dimensional 
HTPEMFC  models  have  been  presented  by  Peng  et  al.  [19,20]. 
They  showed  that  thermal  management  strongly  affects  the  fuel 
cell  performance  and  discussed  key  optimization  parameters  for 
performance  improvements.  Scott  et  al.  [21]  devised  a  one¬ 
dimensional  model.  The  model  was  able  to  satisfactorily  predict 
the  polarization  curve  and  was  used  to  simulate  the  effects  of  cat¬ 
alyst  loading  and  the  Pt/C-ratio  on  fuel  cell  performance.  Ubong 
et  al.  [22]  developed  a  single-channel  three-dimensional  model  in 
which  the  reaction  layer  was  assumed  to  be  infinitely  thin,  and 
the  electrochemical  reactions  were  described  using  an  agglomer¬ 
ate  approach.  A  complete  three-dimensional  model  was  developed 
and  solved  in  [23],  highlighting  reaction-layer  kinetics.  Shamardina 
et  al.  [24]  presented  a  simple  and  quickly  solvable  steady-state, 
isothermal,  pseudo-two-dimensional  model  that  accounted  for 
crossover  effects.  Another  analytical  FITPEMFC  model,  published 
by  Kulikovsky  et  al.  [25],  discussed  important  basic  kinetic  and 
transport  parameters.  A  two-dimensional  isothermal  model  was 
presented  by  Sousa  et  al.  [26].  They  treated  the  reaction  layer 
as  spherical  catalyst  agglomerates  with  porous  inter-agglomerate 
spaces.  The  model  was  used  to  study  the  influence  of  the  reaction 
layer  properties  on  cell  performance.  A  control-oriented,  one¬ 
dimensional  model  was  presented  in  [27],  addressing  the  transient 
responses  of  FITPEMFC.  Wang  et  al.  [28]  investigated  the  transient 
evolution  of  the  carbon  monoxide  poisoning  effect  of  PBI  mem¬ 
brane  fuel  cells  using  a  one-dimensional  model.  Another  work 
that  deals  with  CO  poisoning  and  its  dynamics  is  presented  in 
[29].  Different  FITPEMFC  models  were  presented  in  [30-32],  mainly 
addressed  solid-  and  fluid-(gas)-phase  temperature  and  PBI/FI3PO4 
sol-gel  membrane  behavior.  Structural  mechanics,  namely  local¬ 
ized  fluid-structural  interactions  (LSFI),  were  also  analyzed  [33]. 

The  current  state-of-the-art  FITPEMFC  technology  offers  signif¬ 
icant  optimization  potential.  Besides  material  demands,  operating 
schemes  are  of  great  importance  in  ensuring  effective  and  effi¬ 
cient  fuel  cell  operation,  decreasing  degradation,  and  reaching 
long  operating  times.  FITPEMFC  modeling  and  simulation  may  help 
achieve  these  goals  faster.  This  work  presents  a  general-purpose, 
large-scale  FITPEMFC  model  focusing  on  the  behavior  of  various 
quantities  under  given  operating  conditions.  The  simulated  cur¬ 
rent  density  and  fluid-flow  distribution  should  help  to  optimize 
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Table  1 

Spatial  dimensions  (HTPEMFC  assembly/computational  subdomains)  and  MEA  data  [34]. 
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Symbol 

Value 

Unit 

Note 

Reference 

leu 

0.001 

m 

Thickness 

Ibpp 

0.005 

m 

Thickness 

Ig  dl 

400 x 10“6 

m 

Thickness  (uncompressed) 

Irl,o 

30 x 10-6 

m 

Thickness  (anode  side) 

Irl,c 

40  x 10-6 

m 

Thickness  (cathode  side) 

l  MEM 

152  x 10-6 

m 

Thickness  (uncompressed) 

Qmea 

0.005 

m2 

MEA  active  surface  area 

X 

32 

PBI/H3PO4  doping  level 

[8] 

mPt,a 

0.01 

kg  m-2 

Anode-side  platinum  loading 

[9,10] 

mPt,c 

0.0075 

kg  irr2 

Cathode-side  platinum  loading 

[9,10] 

/pt/C,a 

0.3 

Anode-side  Pt/C  ratio 

/Pt/C,c 

0.3 

Cathode-side  Pt/C  ratio 

the  gas  flow  channel  design  for  HTPEMFC  applications.  As  for  the 
energy  transport,  different  source  term  contributions  are  calcu¬ 
lated.  Results  are  compared  to  experimental  tests  to  make  this 
work  complete.  Solver  settings,  soft-,  and  hardware  requirements 
are  reported  in  great  detail. 

The  work  in  [34]  is  summarized  below.  Theoretical  deriva¬ 
tions  and  experimental  investigations  were  performed  using  a 
HTPEMFC.  Further,  the  electrochemical  behavior  was  highlighted. 
The  maximum  theoretical  thermodynamic  potential  and  maximum 
theoretical  cell  efficiency  was  calculated  using  thermodynamics 
and  compared  to  the  values  of  a  LTPEMFC.  Next,  a  complete  cell 
assembly  (same  as  in  Fig.  2)  was  analyzed  using  EIS  to  deduce 
assembly  resistance  contributions,  followed  by  voltage  loss  esti¬ 
mations.  The  last  section  elucidated  conductivity  measurements 
and  theoretical  calculations  on  the  PBI/H3PO4  sol-gel  membrane 
composition  (BASF  Celtec®-P  Series  2000  MEA).  Some  modeling 
parameters  pertinent  to  HTPEMFCs  are  directly  taken  from  [34]  and 
serve  as  a  basis  for  modeling  and  simulation  (see  Tables  1-3). 

2.  Modeling  aspects 

2.1.  Model  geometry 

The  three-dimensional  model  geometry  used  for  all  computa¬ 
tions  is  given  in  Fig.  1.  It  is  an  exact  representation  of  the  HTPEMFC 
assembly  that  was  used  for  all  experimental  tests,  as  described 
in  [34],  and  includes  gold-plated  copper  current  collectors  (Cu) 
and  four  high-temperature  stable  bipolar  plates  (BPP).  The  gas 
flow  channels  are  machined  into  the  graphite  (six  channel  par¬ 
allel  serpentine  flow-field  with  a  channel-to-land  ratio  of  1. 0/1.0 
(1  x  10-3  m/1  x  10-3  m)).The  gas  flow  channel  depth  is  1  x  10-3  m 
(z-coordinate).  The  MEA  is  sandwiched  between  the  BPPs  and 
includes  the  gas  diffusion  layers  (GDL),  the  reaction  layers  (RL),  and 
a  PBI/H3PO4  sol-gel  membrane  (BASF  Celtec®-P  Series  2000  MEA). 
The  exact  spatial  dimensions  of  the  components  can  be  found  in 
Table  1. 

2.2.  Transport  equations 
22 A.  Momentum  transport 

The  continuity  equation  and  the  incompressible  Navier-Stokes 
equations  (density  assumed  to  be  constant  or  nearly  constant)  are 
solved  to  account  for  the  laminar  gas  flow  and  pressure  distribu¬ 
tion  within  the  anode-  and  cathode-side  gas  flow  channels.  Eq.  (1) 
accounts  for  advection  momentum  flux  and  momentum  imparted 
due  to  pressure  and  viscosity. 

v-u=° 

p  •  (u  •  V)  ■  u  =  V  ■  (-P  ■  I  +  r]  ■  (Vu  +  (Vu)r))  ^  ‘ 

The  Brinkman  equations  (2)  describe  the  gas  flow  within  porous 
media,  i.e.,  GDL  and  RL.  This  mathematical  model  extends  Darcy’s 


law  to  include  a  term  that  accounts  for  the  viscous  transport  in 
the  momentum  balance,  and  it  treats  both  the  pressure  and  the 
flow-velocity  vector  as  independent  variables. 


V  u  —  — 
P 


(J+Sm)  ■  u  =  V  -  (-P./+1.  (^.(Vu  +  fVu)7')-  (|  fj)  (V  u)  /)) 


(2) 


In  Eqs.  (1 )  and  (2),  u  is  the  velocity  vector,  and  P  the  pressure. 


2.2.2.  Mass  (species)  transport 

The  gas  flow  is  predominantly  convective  in  the  gas  flow  chan¬ 
nels,  whereas  it  is  predominantly  diffusive  within  the  porous 
media.  Mass  (species)  transport  is  solved  using  Stefan-Maxwell  dif¬ 
fusion  in  the  convection  application  mode  accounting  for  hydrogen 
(H2 )  and  water  (H20)  at  the  anode  and  oxygen  (02 ),  H20,  and  nitro¬ 
gen  (N2 )  at  the  cathode  (3).  The  reaction  rates  appear  as  source/sink 
terms  S(JL>-  Note  that  the  amount  of  water  in  both  gas  streams  at  the 
inlets  is  small.  This  is  implied  by  our  laboratory  hardware,  but  it 
should  not  be  neglected. 


V  •  p  •  coi  •  u  -  p  •  coi 


J=1 


— 


0) 


2.2.3.  Charge  transport 

Two  Poisson  equations  are  used  to  evaluate  charge  transport  (4). 
The  solid-phase  potential  (subscripts)  distribution  is  solved  within 
the  gold-plated  copper  current  collectors,  the  BPPs,  GDL,  and  the 
RL.  The  membrane-phase  potential  (subscript  m)  is  solved  within 
the  RL  in  the  membrane.  Both  equations  are  coupled  through  their 
current  source/sink  terms. 

— V  •  (crs  •  V0S)  =  —  5^  ... 

-V  •  (crm  •  V0m)  =  +5^  1  ; 


2.2.4.  Energy  transport 

The  thermal  behavior  is  described  using  a  two-equation  system 
(5),  separately  accounting  for  the  solid-phase  (Ts)  and  the  fluid- 
(gas)-phase  temperatures  (T/-),  as  introduced  by  our  group  in  2007 
for  HTPEMFC  models  [30,31  ].  This  is  necessary  because  of  the  large 
temperature  difference  between  the  feed-gas  temperature  (fluid- 
(gas)-phase)  and  the  cell-operating  temperature  (solid-phase)  and 
represents  the  thermal  interactions  between  the  two  phases.  This 
behavior  was  confirmed  by  performing  segmented  temperature 
measurements.  The  solid-phase  temperature  distribution  is  calcu¬ 
lated  for  the  gold-plated  copper  current  collectors,  the  BPPs,  the 
solid  matrix  of  the  porous  media  and  the  membrane.  The  fluid- 
(gas)-phase  temperature  distribution  is  solved  within  the  gas  flow 
channel  and  within  the  porous  media.  Both  equations  are  cou¬ 
pled  through  their  source/sink  terms  and  account  for  possible  heat 
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Fig.  1.  (a)  CAD  representation  of  the  HTPEMFC  assembly;  (b)  detailed  view  of  the  MEA  and  electrochemical  reactions. 


transfer  between  the  solid-  and  the  fluid-(gas)-phase  using  a  volu¬ 
metric  (intrinsic)  heat  transfer  coefficient  ( hTs/Tf ).  In  general,  this 
coefficient  is  a  function  of  the  morphology  of  the  porous  media. 
As  stated  in  [35],  typical  values  for  metal  foams  vary  from  2  x  104 
to  2  x  105Wm-3  I<-1  for  a  porosity  between  0.7  and  0.95.  In  this 
work,  a  similar  value  is  used. 

V  •  (—k  ■  VTS)  =  STs 

V  •  (-k  •  VTf)  =  STf  —  p  Cp  u  ■  V7j  lDj 


2.3.  Source/sink  term  couplings 


Mass  transport  source/sink  terms  are  non-zero  only  in  the  RL 
subdomains.  The  reaction  rates  for  H2,  02,  and  H20  are  described 
with  Eq.  (6). 


S(jh  ~  < 


~ Ja 


Jc 


-Jc- 


Mh2 

Ko 

2  F 


(6) 


Agglomerate  structure  models  have  been  developed  to  describe 
LTPEMFC  and  PAFC  behavior,  among  others  [36].  In  this  work,  an 
agglomerate  model  is  used  to  simulate  electrochemical  half-cell 
reactions.  The  LTPEMFC  model  from  COMSOL  Multiphysics®  [37] 
and  the  equations  presented  in  [38]  serve  as  a  foundation  and  are 
adopted  for  HTPEMFC  modeling  and  simulation  by  accounting  for 
gas  diffusivity  and  gas  solubility  in  H3P04,  and  the  amorphous 
phase  of  H3PO4,  respectively.  It  is  referred  to  [39]  for  a  detailed 
description  of  the  analytical  solution  of  a  diffusion  reaction  prob¬ 
lem  in  a  spherical  porous  particle.  In  this  work,  the  reactions  are 
assumed  to  occur  in  electrolyte-filled  agglomerate  zones.  These 
zones  contain  channels  of  the  amorphous  phase  H3PO4  into  which 
the  gases  dissolve.  Gas-transport  limitations  of  diffusion  processes 
in  the  agglomerates  are  characterized  by  an  effectiveness  factor 
and  the  Thiele  modulus  (dimensionless  group).  The  characteristic 
length  scale  is  the  radius  of  the  agglomerate  ragg.  Eq.  (7)  describes 
the  volumetric  current  density  at  the  anode. 


The  transfer  coefficients  aa  and  ac  in  Eqs.  (7)  and  (8)  depend 
on  many  parameters,  including  temperature.  In  this  work,  at  the 
anode-side,  a  well-recognized  value  of  0.5  was  assumed,  see  e.g. 
[26].  At  the  cathode-side,  the  reported  values  vary  within  a  certain 
range.  In  [26],  a  value  of  0.73  best  described  the  behavior  of  the  base 
system.  Huang  et  al.  [41]  reported  a  value  of  0.67  at  150°C  for  02 
reduction  on  Pt  in  85%  H3PO4.  Kulikovsky  et  al.  [25]  discussed  the 
Tafel  slope  and  used  a  value  of  0.777  for  the  transfer  coefficient.  In 
this  work,  a  cathode-side  transfer  coefficient  of  0.89  returned  the 
best  results. 

As  for  solid-phase  energy  transport,  the  source/sink  terms  are 
given  in  Eq.  (9)  and  account  for  ohmic  and/or  protonic  heat¬ 
ing  (membrane  heating),  irreversible  reaction  heat,  and  reaction 
entropy  for  the  following  subdomains:  the  gold-plated  copper  cur¬ 
rent  collector  and  bipolar  plates,  GDL,  anode-side  RL,  cathode-side 
RL,  and  the  membrane  (top  to  bottom). 


STs  =  < 


Ts/Tf  •  (Ts  -  Tf ) 


1 

—  -  h 

a4  .2 

-Jr  +  +  }a  •  Va  -  hTs/Tf  ■  ( Ts  -  Tf ) 

<7j  <7m  1 

h  ,  ,  j  ,  j  ASC  •  Ts  ,  / t*  t1  ^ 

^  ^  Jc'k  +]c ' ~^rr ~ hr>'Tf ' (Ts ~  f] 

!k 

(7m 


(9) 


Eq.  (10)  introduces  the  heat  source/sink  term  for  the  fluid-(gas)- 
phase  temperature  within  the  GDL  and  RL  subdomains. 


$rf  =  hTs/Tf  •  (Ts  -  Tf) 


(10) 


2.4.  PBI/H3PO4  sol-gel  membrane  modeling 

The  following  equations  describe  the  behavior  of  a  highly  doped 
PBI/H3PO4  system  of  approximately  30-35  mol  of  H3P04  per  PBI 
repeat  unit.  As  mentioned  in  [34],  the  properties  of  the  membrane 
are  significantly  affected  by  the  polymer  molecular  structure  and 


ja  =  -(1  -  £rL)  •  3  •  2  •  F  ■  Dh2 ,pBI/(X-2)H3P04/ragg  '  (CH2,PBI/(X-2)H3P04  -  cu2  *  '^a)  •  (1  -  I<a  •  coth(/Ca)) 

Ka  =  ^/'pBI/(X-2)H3P04,a  '  Qa  •  1/2  ■  F  ■  •  Dh2,PBI/(X-2)H3P04  '  ragg 


The  cathode  volumetric  current  density  is  described  using  Eq.  (8). 

jc  =  (1  -  £rL)  •  3  •  4  •  F  ■  Do2iPBl/(X-2)H3P04/ragg  *  c02,PBI/(X-2)H3P04  •  (1  -Kc-  coth (Kc)) 

Kc  =  a/7PBI / (X-2 )H3 P04 , c  •  ac  •  1/4  •  F  •  cj2  •  Do2ipbi/(X-2)H3P04  •  e  ac  F/R  Tf  Vc  •  ragg 
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the  preparation  process  itself.  Because  two  H3P04  molecules  are 
bonded  to  PBI,  X-2  molecules  remain  free  and  tend  to  form  an 
amorphous  phase  within  the  PBI/H3PO4  sol-gel  membrane  [7]. 
It  is  expected  that  the  final  membrane  initially  consists  of  more 
than  85  wt.%  H3P04.  The  amorphous  phase  H3P04  volume  fraction 
£pbi/(x-2)h3po4  within  the  membrane  is  calculated  using  Eq.  (11). 


£PBI/(X-2)H3P04  = 


(  Mpbi/Mh3po4  + 

V  ^2  J 


(11) 


It  is  accepted  that  the  amorphous  phase  H3P04  volume  fraction 
contributes  to  the  high  membrane  conductivity  via  a  Grotthuss  pro¬ 
ton  switching  mechanism.  The  strong  temperature  dependency  of 
the  membrane  conductivity  can  be  described  using  an  Arrhenius 
approach  [7]. 

a  =  a0(khX)/Ts  •  e(-A£a^*)/^)  (12) 

The  pre-exponential  term  a0  is  assumed  to  be  independent  of 
the  cell  operating  temperature  and  decreases  with  higher  doping 
levels  X.  The  concentration  of  the  mobile  species  in  Eq.  (13)  should 
change  with  the  doping  level  [7]. 

^  •  a  ■  v0  •  d2  •  c  •  eAS+Asf/R  (13) 


The  activation  energy  of  the  conductivity  depends  on  multiple 
factors,  including  the  membrane  doping  level  and  the  polymer 
backbone  structure.  Only  minor  amounts  of  valuable  data  are  avail¬ 
able  in  the  literature  for  PBI/H3P04  sol-gel  membranes  [11,12].  In 
this  work,  the  obtained  values  for  the  activation  energy  and  pre¬ 
exponential  factor  are  taken  from  [34]. 


2.5.  Additional  modeling  equations 


In  [40],  it  was  demonstrated  that  the  logarithmic  exchange  cur¬ 
rent  density  of  oxygen  reduction  at  platinum-interfaced  PBI/H3P04 
linearly  increases  with  the  amorphous  phase  H3P04  volume  frac¬ 
tion.  From  this  data,  the  authors  of  that  report  estimated  the 
exchange  current  density  for  concentrated  H3P04  and  reported 
consistent  values.  In  this  work,  Eq.  (14)  is  used  to  calculate  the 
cathode-side  exchange  current  density  ipBI/(X_2)h3po4  c* 

1PBI/(X-2)H3P04,c  =  IH3P04,c  ‘  }  q4.  16  (  i-ePBI/(X-2)H3P04 ) 

In  [41],  it  was  found  that  the  linearity  of  the  logarithmic 
current  density  and  the  reciprocal  temperature  plot  indicates 
Arrhenius  behavior  with  an  apparent  activation  energy  of 
41 840 J  mol-1  (85  wt.%  H3P04).  Reported  values  are  in  the  range  of 
0.01-0.02  Am-2  for  typical  operating  temperatures.  In  the  present 
work,  the  exchange  current  density  for  concentrated  H3  P04  i9 .  p04  c 
is  calculated  with  Eq.  ( 1 5 ).  A  factor  of  1  x  1 04  is  used  to  convert  into 
SI  units  (Am-2). 

i°3p04?c  =  1  X  104  X  io^-0-491-2193  1/7/)  (15) 

For  the  anode  half-cell  reaction,  the  exchange  current  density  is 
taken  to  be  1  x  108  times  the  cathode-side  exchange  current  den¬ 
sity,  a  value  that  is  consistent  with  many  other  published  works, 
see  e.g.  [18].  The  effective  surface  area  of  the  RL  (16)  is  a  function  of 
the  catalyst  surface  area  per  unit  mass  of  the  catalyst  particle,  the 
platinum  loading,  and  the  thickness  of  the  RL  [42]. 


ai  =  aPt,i  • 


mPt,i 

Irlj 


(16) 


Gallart  [43]  summarized  the  available  catalyst  surface  area  per  unit 
mass  of  the  catalyst  particle  data  and  coupled  it  to  the  platinum- 


Table  2 

Modeling  constants  and  parameters. 


Symbol 

Value 

Unit 

Reference 

A  Ea 

18,484 

J  mol-1 

[34] 

F 

96,485 

As  mol-1 

^BPP/CH  ,a 

15 

W  nrr2 K-1 

^BPP/CH,c 

40 

W-2  K-1 

hTs/Tf 

1  x  105 

W  nrr3 K_1 

[41] 

k 

3.16  x  10“8 

kcu 

400 

Witt1  K"1 

[37] 

l<  BPP 

20 

Wm-1  K"1 

/<GDL,x,y 

2.5 

Wm-1  K"1 

/<GDL,z 

0.25 

Wm-1  K"1 

fcp.GDL 

5  x  10-13 

m2 

l<p,  RL 

1  x  10-13 

m2 

1<PBI 

0.45 

Wm-1  K"1 

1<RL 

0.3 

Wm-1  K"1 

Mh2 

0.002 

kg  mol-1 

Mq2 

0.032 

kg  mol-1 

Mn2 

0.028 

kg  mol-1 

Mh2o 

0.018 

kg  mol-1 

Mh3po4 

0.098 

kg  mol-1 

Mpbi 

0.308 

kg  mol-1 

ragg 

1  x  10-9 

m 

[43] 

R 

8.314 

J  mol-1 1<-1 

aa 

0.5 

OLc 

0.89 

g0 

&GDL 

0.78  (uncompressed) 

ln( cr0  ■  Ts) 

9.5211 

SKcm-1 

[34] 

&CU 

1  x  104 

Sm_1 

[37] 

ctbpp 

3600 

Sm4 

G  GDL  ,x,y,z 

220 

Srrr1 

[40] 

gs,rl 

450 

Sm_1 

[38] 

Gm,  RL 

13 

Sm4 

[38] 

VH2 

7.07  x  10“6 

m3  mol-1 

v02 

16.6  x  10“6 

m3  mol-1 

Vn2 

17.9  x  10“6 

m3  mol-1 

Vh2o 

12.7  x  10“6 

m3  mol-1 

Pc 

2000 

kg  m-3 

[37] 

PpBI/(X-2)H3P04 

1698 

kg  m-3 

[40] 

Ppt 

21,500 

kg  m-3 

[37] 

to-carbon  ratio,  as  shown  in  Eq.  (17). 

apt,,  =  7.401  x  105  /p4t/c  ,-1.811  x  106  /P3t/C  ,-+1.545  x  106  /p2t/c  j 

-6.453  x  105  •/Pt/c,i+2.054  x  105  (17) 

The  local  RL  overpotential  r)i  is  defined  as  the  difference  between 
the  solid-  and  membrane-phase  potentials.  This  potential  differ¬ 
ence  drives  the  cell  current,  keeping  the  electrochemical  half-cell 
reactions  continuous.  It  is  calculated  with  Eqs.  (18)  and  (19). 

FJa  —  0s  —  0m  0  (18) 

dc  =  (Ps  0m  —  (19) 

The  maximum  equilibrium  potential  with  respect  to  temperature 
and  partial  pressures  is  calculated  using  the  Nernst  equation,  as 
discussed  in  [34]  for  a  HTPEMFC. 

Accounting  for  gas  diffusivity  and  solubility  in  the  electrolyte  is 
crucial  for  HTPEMFC  modeling,  as  mentioned  in  the  few  currently 
available  studies  [  1 8,26,40].  In  [40],  the  authors  found  that  the  oxy¬ 
gen  diffusivity  increases  with  increasing  PBI/(X-  2)H3P04  and  that 
it  does  not  significantly  vary  from  the  data  available  for  hot  con¬ 
centrated  H3P04.  In  the  present  work,  a  similar  expression  is  used, 
as  presented  in  [18].  A  Bruggeman  relation  with  an  exponent  of  1.8 
is  used  to  account  for  the  agglomerate  structure  (20)  [40]. 

Ho2,PBI/(X-2)H3P04  Ho2,H3P04  ’  (£pBI/(X-2)H3P04  '  ^Rl)  (20) 

The  oxygen  solubility  is  generally  higher  (about  four  times  higher 
than  expected  for  pure  H3P04  under  given  conditions  [14])  than 
the  values  reported  for  hot  concentrated  H3P04.  In  [40],  the  authors 
stated  that  these  higher  values  must  be  related  to  the  presence  of 
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PBI.  Cheddie  and  Munroe  [18]  related  the  values  of  hot  concentrated 
H3PO4  to  the  values  presented  by  Liu  et  al.  [40]  using  Eq.  (21 ). 

^02,PBI/(X-2)H3P04  =  ePB^/(* * * * 5X-2)H3P04  '  ((^02,H3P04 ) 

+5.79  •  (1  -  £Jbi/(x-2)h3po4)) 

Data  for  the  hydrogen  oxidation  reaction  at  the  anode-side  are 
not  readily  available,  so  the  hydrogen  diffusivity  is  taken  to  be 
two  times  the  oxygen  diffusivity  and  the  hydrogen  solubility 
is  taken  to  be  approximately  four  times  the  oxygen  solubility. 
The  same  behavior  is  assumed  as  in  a  water  system  [18].  On 
the  other  hand,  one  has  to  note  that  Li  et  al.  [14]  reported 
solubility  values  of  1.6  x  10-5  mol  cm-3  bar-1  for  hydrogen  and 
1.9  x  10-5  mol  cm-3  bar-1  for  oxygen  in  PBI  membranes. 


£)H2,PBI/(X-2)H3P04  =  2  •  £)02,PBI/(X-2)H3P04  (22) 

^h2,pbi/(x-2)h3po4  =  4.44  •  Hq2 ,pbi/(x-2)h3po4  (23) 

The  gas  diffusion  coefficient  and  gas  solubility  in  the  amorphous 
phase  H3PO4  must  also  be  related  to  the  reported  values  for  hot 
concentrated  H3P04.  Klinedinst  et  al.  [44]  highlighted  that  oxy¬ 
gen  diffusivity  and  solubility  in  hot  concentrated  H3PO4  exhibit 
exponential  reciprocal  temperature  dependencies  over  sufficiently 
small  temperature  ranges.  The  activation  energy  for  oxygen  dif¬ 
fusion  and  the  enthalpy  of  the  solution  vary  with  concentration. 
In  this  work,  the  transport  properties  of  oxygen  in  concentrated 
H3PO4  are  related  to  the  temperature  and  the  acid  concentration, 
as  presented  in  [18].  Similar  empirical  equations  were  used  in  [26]. 


_ q  ((-192.55.m2 

^02,H3P04  =  1x10  e  3 


H,PO4+323-55mH3PO4-125.61)+62010.m^pO4 


Table  3 

Base  case  operating  parameters  [34]. 


T0 o  (ambient)  (°C) 

21 

Tceu  (solid-phase)  (°C) 

160 

Anode  side 

Cathode  side 

Gas 

Hydrogen 

Air 

Flow  rate  (lmin-1) 

Stoichiometry  1.35 

Stoichiometry  2.5 

Gas  inlet  temperature  (°C) 

21 

21 

Gas  humidity  (%  rh) 

0.2 

2 

Outlet  pressure  (Pa) 

1.01325  x  105 

1.01325  x  105 

2.6.  Solid-  and  fluid- (gas) -phase  properties  and  material 
correlations 


A  woven-type  GDL  (similar  to  E-tek  ELAT®  products  [45])  is 
considered  herein.  It  consists  of  a  void  volume  fraction  (porosity; 
superscript  0)  and  a  solid-phase  volume  fraction  (subscript  s),  as 
indicated  by  Eq.  (30). 

£gdl  +  sgdl  =  1  (30) 

The  RL  structure  is  more  complex.  Different  volume  fractions  are 
needed  to  calculate  the  effective  RL  properties,  e.g.,  thermal,  elec¬ 
trical  or  protonic  conductivities.  The  volume  fraction  of  platinum 
and  carbon  (subscript  Pt/C)  is  calculated  with  Eq.  (31),  taking  its 
thickness  and  the  catalyst  properties  into  account  [46]. 


£pt/c  = 


raPt 

Irl  /pt/c 


1  -fpt/c\ 
Pc  ) 


(31) 


-1 05503  •i7iH3po4+40929/7y) 


(24) 


Ho2,h3po4  =  1x10  1  •  e 


((257.13.m2  -431.08.mH3Po4+178.45)+-93500.m2  +156646.mH3Po4-64288/r/) 


(25) 


The  H3PO4  concentration  within  the  MEA  is  expected  to  change 
during  cell  operation  (compare  with  Fig.  13).  In  [42],  a  correlating 
equation  for  H3PO4  vapor  pressure  is  given  (80-101  wt.%  acid  in 
the  range  of  1 30-1 70  °C).  From  their  experimental  data,  Souza  et  al. 
[26]  generated  an  equation  that  coupled  concentration  and  water 
vapor  partial  pressure  (26). 


*h3po4  = 


ln(xH2o  •  P)  +  2765.1/7/  -  22.002 
-4121.9/7) +  2.5929 


(26) 


A  similar  equation  was  used  by  Choudhury  et  al.  [36]  to  corre¬ 
late  the  H3PO4  concentration  and  equilibrium  vapor  pressure  at 
a  given  temperature  (equilibrium  between  the  humidity  and  the 
H3PO4/H2O  solution).  The  mole  fraction  of  H3PO4  is  converted  into 
mass  fraction  with  Eq.  (27),  similar  to  [42]. 


In  this  model,  the  fraction  of  the  volume  occupied  by  the  amor¬ 
phous  phase  H3PO4  inside  the  agglomerate  is  calculated  with  Eq. 
(32)  [46]. 

_  rapt  /  /pbi/(x-2)h3po4  \ 

^  4  kl-fpt/C  y  (1  -/pbI/(X-2)H3P04)  •  Ppbi/(x-2)h3po4  J 

(32) 

Finally,  the  RL  porosity  can  be  calculated  with  Eq.  (33). 

SRL  =  1  -  £Pt/C  -  £pBI/(X-2)H3P04  (33) 

The  binary  diffusion  coefficients  Dy  are  calculated  for  all  pairs  of 
species  in  the  fluid-(gas)-phase  mixture  using  Eq.  (34)  and  are 
known  to  vary  with  pressure  and  temperature. 


raH3po4 


136  Xh3po4 
111  -Xh3po4  +  25 


(27) 


Finally,  the  hydrogen  and  oxygen  agglomerate  concentration  is 
calculated  using  Henry’s  law  (28)  and  (29),  incorporating  the 
partial  pressures  of  hydrogen  and  oxygen  (gas/electrolyte  (amor¬ 
phous  phase  H3PO4)  interface).  Note  that  different  solubility  units 
are  used  in  the  literature,  e.g.,  mol  m-3  atm-1,  mol  m-3  Pa-1  or 
molm-3  bar-1. 


cH2,PBI/(X-2)H3P04  = 

Px  h2 

(28) 

Wh2,PBI/(X-2)H3P04 

c02,PBI/(X-2)H3P04  = 

P-Xo2 

(29) 

^02,PBI/(X-2)H3P04 

Dij  =  k. 


M+++3) 


1  1 

Wi  +  Wj 


(34) 


The  effective  porous  media  binary  diffusivities  Df  in  Eq.  (3)  are 
calculated  using  the  respective  void  volume  fraction  and  a  Brugge- 
man  relationship  with  an  exponent  of  1.5,  which  is  consistent  with 
other  works. 


heff  —  n..  .  ( F0 
i7  “  UV 


+  .5 


GDL.RL' 


(35) 


The  densities  of  the  fluid-(gas)-phase  mixture  within  the  gas  flow 
channel  and  the  porous  media  are  calculated  using  the  mole  frac¬ 
tion  and  the  molar  mass  of  the  gas  species. 


P 

RTf 


(36) 
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The  fluid-(gas)-phase  properties,  namely  the  density,  thermal 
conductivity,  heat  capacity  and  dynamic  viscosity,  are  taken  from 
the  material  library  provided  by  the  software  used  and  depend  on 
the  fluid-(gas)-phase  temperature. 

The  thermal  conductivity  of  the  membrane  /<mem  is  calculated 
using  the  thermal  conductivities  of  PBI  and  amorphous  phase 
H3PO4  (amorphous  phase  H3PO4  treated  as  concentrated  H3PO4). 

^MEM  =  f<PBI  •  (1  -  £PBI/(X-2)H3P04)  +  ^PBI/(X-2)H3P04  *  £PBI/(X-2)H3P04 

(37) 

Based  on  the  data  provided  by  Turnbull  [47],  the  thermal  con¬ 
ductivity  data  is  extrapolated  for  higher  temperatures  and  H3PO4 
concentrations  with  Eq.  (38). 


assumed  to  be  85wt.%  at  room  temperature  and  pressure 
(RTP). 

(8)  The  membrane  is  considered  to  be  a  system  of  PBI,  H3PO4 
and  H20  only  (molecular  H3P04,  other  H3PO4  forms  are 
neglected). 

(9)  Low  velocities  (low  Reynolds  numbers)  are  expected  in  the 
gas  streams.  Laminar  inflow  conditions  can  be  assumed. 

(10)  All  agglomerates  are  assumed  to  be  geometrically  identical, 
spherical  in  shape,  and  of  the  same  radii.  An  additional  H3P04 
(or  water)  film  surrounding  the  agglomerate  is  not  considered. 

(11)  There  is  no  heat  transfer  towards  the  surroundings. 

(12)  There  is  no  (partial)  H3PO4  flooding  of  the  GDL  and  no  H3PO4 
catalyst  absorption 


'<h3po4  = 


11.727  +  0.01864  •  Ts  -  0.02169  •  cH3po4  -  0.0000338  •  cH3po4  •  Ts 

MO* * 3 4 


•418.68 


(38) 


kpBI/(X-2)H3P04  ~ 

The  remaining  gas  mixture  properties  were  calculated  using 
average-based  mole  or  mass  fractions.  Other  modeling  parameters 
can  be  found  in  Table  2.  Some  of  these  values  (or  their  orders  of 
magnitude)  were  provided  by  their  respective  manufacturers. 

2.7.  Boundary  conditions 

Boundary  conditions  are  crucial  for  accurate  simulation  results. 
They  are  precisely  applied  according  to  the  experimental  setup 
(Table  3).  The  gas  composition  (mass  fractions),  gas  temperature, 
and  gas  velocity  (mass  flux)  are  chosen  at  the  gas  flow  channel 
inlets.  The  inlet  flow  should  be  fully  developed  when  entering  the 
cell.  An  additional  pressure  variable  (weak  contribution/additional 
degrees  of  freedom  (DOF))  is  used  to  define  a  laminar  inflow 
condition.  Reference  concentrations  are  deduced  from  the  inlet 
conditions.  At  the  outlets,  pressure  and  convective  flux  bound¬ 
ary  conditions  are  used.  To  conserve  charge,  the  cell  voltage  is 
given  and  the  current  density  calculated.  The  cell  operating  voltage 
is  defined  at  the  cathode-side  gold-plated  copper  current  col¬ 
lector.  At  the  anode-side,  the  potential  is  fixed  at  zero.  The  cell 
operating  temperature  (solid-phase)  is  defined  at  the  anode-  and 
cathode-side  gold-plated  copper  current  collector  boundaries.  At 
the  bipolar  plate  walls  (gas  flow  channel  walls),  a  no-slip  boundary 
condition  is  defined  for  momentum  transport,  and  a  heat  transfer 
coefficient  is  used  to  account  for  the  heat  by  the  absorbed  gases 
in  both  energy  transport  equations.  At  all  other  internal  bound¬ 
aries,  symmetry,  continuity  or  insulation  boundary  conditions  are 
used. 

2.8.  Assumptions  and  simplifications 

The  following  assumptions  and  simplifications  are  used: 

(1)  Steady-state  operating  conditions. 

(2)  Continuity  is  prescribed  at  internal  boundaries  (interfaces);  all 
contact  resistances  are  neglected. 

(3)  Some  material  properties  are  assumed  to  be  isotropic  and 
macrohomogeneous,  whereas  others  account  for  the  different 
material  properties  in  the  x-,  y-,  and  z-directions. 

(4)  A  microporous  layer  and  its  influence  on  the  quantities’  behav¬ 
ior  and  distribution  is  not  explicitly  considered  in  this  model. 

(5)  All  product  water  is  assumed  to  be  vaporous  because  of  the 
high  operating  temperature  (no  phase  change)  and  to  leave 
the  cell  assembly  in  vapor  form. 

(6)  Gas  and  water  crossover  through  (or  water  uptake  by)  the 
PBI/H3PO4  sol-gel  membrane  is  neglected. 

(7)  The  initial  concentration  of  the  amorphous  phase  H3PO4 
volume  fraction  inside  the  RL  and  inside  the  membrane  is 


The  listed  assumptions  are  generally  adopted  when  modeling  a 
HTPEMFC  (e.g.,  assumptions  2, 5, 9).  Some  of  the  simplifications  are 
directly  related  to  steady-state  operating  conditions  (assumption 
1 ),  the  modeling  level,  the  purpose,  and  the  outcome  of  the  study 
(e.g.,  contact  resistances  are  mostly  neglected  in  similar  studies 
but  should  be  included  when  modeling  cell  compression,  structural 
mechanics,  and/or  cell  ageing).  Other  simplifications  base  on  mate¬ 
rial  properties  and/or  material  casting  methods  (e.g.,  assumption 
7).  There  is  no  doubt  that  assumptions  4, 6, 1 0,  and  1 2  are  the  major 
drawbacks  of  this  model  setup  and  should  be  updated  in  future 
works.  Assumption  6  is  used  since  the  model  deals  with  steady- 
state  operating  conditions,  high  operating  temperatures  and  low 
vapor  pressure.  Start-up  and  shut-down  periods  are  not  investi¬ 
gated.  Assumption  11  is  discussed  in  more  detail  in  Section  5.4. 

Based  on  the  simulation  results  and  the  observations  made 
during  experimental  test,  it  can  fairly  be  concluded  that  all  assump¬ 
tions  are  acceptable  when  focusing  complete  3D  cell  level  modeling 
and  simulation. 

3.  HTPEMFC  operation 

The  reference  HTPEMFC  in  Fig.  2  includes  a  commercially  avail¬ 
able  PBI/H3PO4  sol-gel  membrane  (Celtec®-P  Series  2000  MEA) 
that  is  stable  at  high  temperatures.  The  available  MEA  data  can 
be  found  in  Table  1.  The  HTPEMFC  assembly  was  temperature- 
controlled  using  heating  elements.  After  steady  state  operating 
conditions  were  achieved,  a  defined  amount  of  hydrogen  and  air 
was  fed  in  a  counter-flow  configuration  using  mass  flow  controllers. 
The  solid-  and  fluid-(gas)-phase  temperature  was  measured  close 
to  the  inlets  and  outlets.  Additionally,  the  pressure  and  relative 
humidity  was  monitored  at  the  in-  and  outlets. 

Hydrogen  at  the  anode-side  and  air  at  the  cathode-side  pass 
the  gas  flow  channels  (flow-field)  and  partially  diffuse  through  the 
porous  media  towards  the  RL,  in  which  the  electrochemical  half¬ 
cell  reactions  take  place  according  to  the  equations  given  in  Fig.  1 . 
The  resulting  voltage  of  a  single  fuel  cell  is  in  the  order  of  1 V.  Power 
can  be  drawn  using  an  external  load  connected  to  both  gold-plated 
copper  current  collectors.  Moreover,  water  and  heat  is  produced  by 
the  exothermic  reactions. 

4.  Computational  aspects 

4.1.  Software  and  hardware 

The  model  geometry  was  created  using  CAD  software  and 
imported  as  a  neutral  standard  (*.iges)  into  commercially  available 
FEM/CFD-software  COMSOL  Multiphysics®  v3.5a  (64  bit  version) 
running  on  a  PC  (double  quadcore  with  64  GB  RAM).  Next,  the  com- 
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Fig.  2.  (a)  Gas  flow  channels;  (b)  HTPEMFC  assembly  used  for  all  experimental  tests. 


plete  3D  geometry  was  discretized  using  355,431  tetrahedral  and 
397,644  prism  elements  (Fig.  3).  Meshing  needed  to  be  performed 
carefully  to  minimize  the  number  of  DOF  and  therefore  the  mem¬ 
ory  requirements.  However,  the  quality  of  the  mesh  was  crucial  for 
accurate  (mesh-independent)  simulation  results.  The  membrane 
subdomain  was  only  meshed  with  tetrahedral  elements  because 
both  bipolar  plates  (and  flow-fields)  are  turned  180°  against  each 
other.  This  results  in  a  different  boundary  mesh  at  both  mem- 
brane/RL  interfaces.  The  mesh  of  the  porous  media  subdomains 
was  refined  to  account  for  the  expected  local  quantity  gradients, 
whereas  a  coarser  mesh  was  used  within  the  gold-plated  copper 
current  collector  and  BPP  subdomains  (only  scalar  variables  had  to 
be  solved  within  these  subdomains).  The  total  number  of  DOF  for 
this  problem  was  12,665,886. 


4.2.  Solution  procedure  and  convergence  behavior 

Solving  a  large  scale  fuel  cell  model  using  finite  elements  is 
not  easy.  A  large  amount  of  memory  is  needed,  even  when  using 
iterative  solvers  (compare  with  Table  4).  Additionally,  one  has  to 
carefully  tune  these  solvers  to  reach  a  converged  solution.  The 
best  possible  initial  conditions  were  generated  using  a  series  of 
parametric  dummy  simulations.  The  problem  was  solved  using 
the  chemical  engineering  and  heat  transfer  modules  provided  by 
the  software.  The  solution  procedure  was  rather  complex  due  to 
the  strong  multiphysical  couplings  of  all  transport  equations,  as 
explained  in  [48].  The  complete  Navier-Stokes  equations  were 
solved  first,  followed  by  the  charge  and  mass  (species)  transport 
equations  and,  finally,  the  solid-  and  fluid-(gas)-phase  tempera- 
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Fig.  3.  (a)  Generated  3D  computational  mesh  inx-y-z-plane;  (b)  mesh  details  iny-z-plane;  (c)  gas  flow  channel  details  (180°  turned  against  each  other). 
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Table  4 

Details  recorded  while  solving  the  HTPEMFC  model  (base  case  operating  conditions/20  A  load  current). 
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Variables  (DOF) 

Solver  settings 

Iterations 

Memory/GB 

Clock  time/s 

Convergence 

uc,  Pc,  PinLchnSc  (2340478) 

PARDISO 

16 

50.4 

19064 

1  x  10“5 

ua,  Pa,  PinLchnSa  (2341345) 

PARDISO 

6 

51.3 

7005 

1  x  10-5 

(pi,  a)i  (4622790) 

Segregated  group  solver/GMRES  and  PARDISO 

7/210 

38.6 

9516 

1  x  10-6 

Ts,  1/(3361273) 

Segregated  group  solver/GMRES  and  PARDISO 

4/36 

17.6 

2731 

1  x  10“4 

ture  distributions  in  an  iterative  fashion.  Table  4  lists  the  details 
recorded  while  solving  the  model  for  base  case  operating  condi¬ 
tions.  It  was  observed  that,  for  such  large  models,  the  system  matrix 
factorization  generally  takes  longer  than  the  solution  process  itself. 
Moreover,  the  number  of  iterations  and  calculation  time  required 
to  meet  the  convergence  criteria  mainly  depended  on  the  cell  oper¬ 
ating  voltages.  The  time  required  for  calculating  a  typical  I-V  curve 
was  roughly  48  h  of  clock  time. 

5.  Results  and  discussion 

5.1.  HTPEMFC  overall  performance 

Fig.  4  displays  the  simulated  HTPEMFC  overall  performance  in 
terms  of  the  polarization  curve  at  the  base  case  operating  condi¬ 
tions  found  in  Table  3.  The  curve  shows  typical  HTPEMFC  behavior. 
The  theoretical  open  circuit  potential  is  much  higher  than  the  mea¬ 
sured  one  (0.168  V  about  the  measured  value),  possibly  due  to  gas 
crossover  and/or  material  imperfections.  At  high  cell  voltages,  a 
sharp  voltage  drop  due  to  activation  losses  (kinetics)  is  observed.  In 
contrast,  the  voltage  drop  is  quite  low  (0.031  V  per  5  A)  in  the  ohmic 
region.  At  a  typical  cell  current  of  20  A  (160°C),  the  cell  voltage 
is  0.5957  V.  At  different  operating  temperatures,  the  cell  returned 
0.5724  V  ( 1 40  °C)  and  0.61 64  V  ( 1 80  °C).  For  a  20  A  load  current,  the 
power  density  is  2382  Wm-2  (160 °C),  2465  Wm-2  (180 °C),  and 
2289  W  m-2  ( 1 40  °C),  respectively.  The  HTPEMFC  model  performed 
better  at  higher  operating  temperatures  and  worse  at  lower  oper¬ 
ating  temperatures.  This  is  consistent  with  the  observations  made 
during  all  experimental  investigations. 

The  different  contributions  to  the  total  HTPEMFC  resistance  are 
discussed  in  [34].  The  total  ohmic  voltage  loss  (from  the  gold-plated 
copper  current  collector  to  RL)  is  calculated  to  be  20.435  mV  for  a 
20  A  load  current.  Under  base  case  operating  conditions  and  a  20  A 
load  current,  the  mean  overpotential  at  the  anode-  and  cathode- 
sides  was  calculated  to  be  0.001  V  and  -0.485  V,  respectively. 


The  developed  model  was  not  able  to  return  the  measured  I-V 
curve  exactly  and  tended  to  slightly  overpredict  the  cell  perfor¬ 
mance  for  low  current  densities,  though  good  agreement  is  found 
for  higher  load  currents.  The  calculated  cathode-side  exchange  cur¬ 
rent  density  (Eqs.  (14)  and  (15))  was  tuned  with  a  constant  in  such 
a  way  that  the  model  returned  approximately  4000  A  nrr 2  at  a  typ¬ 
ical  cell  voltage  of  0.6  V.  This  tuned  value  was  not  changed  further 
for  other  cell  operating  voltages  or  operating  conditions.  This  indi¬ 
cates  that  the  real  cathode-side  exchange  current  density  may  be 
higher  than  the  value  calculated  with  Eqs.  (14)  and  (15). 

5.2.  Velocity  and  pressure  distribution 

Momentum  transport  is  calculated  in  the  gas  flow  channels 
and  within  the  porous  media.  Hydrogen  and  air  (oxygen)  streams 
through  the  gas  flow  channels  and  diffuses  through  the  void  vol¬ 
ume  of  the  GDL  and  reaches  the  RL.  Additionally,  the  gases  have  to 
diffuse  into  the  electrolyte,  as  explained  above.  At  stoichiometric 
flow  rates  for  a  20  A  load  current,  the  mean  gas  flow  channel 
velocity  was  calculated  to  be  2.1584  m  s-1  at  the  cathode  side  and 
0.4914 ms-1  at  the  anode  side.  For  all  operating  conditions  under 
consideration,  the  mean  velocity  is  much  lower  at  the  anode  side 
than  on  the  cathode  side.  Because  of  the  total  flow  rates,  the  mean 
velocity  increases  linearly  when  drawing  more  power  from  the 
cell.  Moreover,  the  velocity  towards  the  outlet  slightly  increases 
as  oxygen  is  consumed  and  water  is  produced,  but  no  significant 
change  in  velocity  was  observed  within  the  anode-side  gas  flow 
channel. 

Within  the  porous  media,  the  mean  velocity  values  were  calcu¬ 
lated  to  be  several  orders  of  magnitude  lower  than  those  in  the  gas 
flow  channels  (e.g.,  0.001 91 6  m  s-1  within  the  cathode-side  porous 
media  at  stoichiometric  flow  rates  for  a  20  A  load  current).  Fig.  5 
depicts  the  local  gas  velocities  within  the  middle  of  the  GDL. 

When  taking  a  closer  look  at  the  velocity  distribution  within 
the  porous  media,  it  is  seen  that  much  higher  velocity  values  are 


Fig.  4.  Measured  and  simulated  HTPEM  fuel  cell  performance  at  base  case  operating 
conditions. 
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Fig.  5.  Local  anode-  and  cathode-side  GDL  velocity  along  the  x-axis  (y/ymax  =  1/2, 
z  =  -150  x  10-6  m)  at  different  stoichiometric  flow  rates  without  drawing  current 
from  the  cell. 
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Fig.  6.  Local  anode  and  cathode-side  GDL  pressure  loss  along  the  x-axis  (y/ym ax  =  1/2, 
z=-150  x  10-6  m)  at  different  stoichiometric  flow  rates  without  drawing  current 
from  the  cell. 


Fig.  8.  Cathode-side  RL  current  density  distribution  in  along  the  x-axis  (plot  A- A', 
z  =  -379  x  10_6m)  at  base  case  operating  conditions  for  5  A,  10  A,  15  A,  20  A,  and 
25  A  load  currents. 


observed  in  the  region  under  each  180°  bend.  This  is  due  to  the 
fact  that  a  pressure  difference  (e.g.,  50  Pa  at  a  15  A  load  current) 
exists  at  the  same  locations,  as  can  be  seen  in  Fig.  6.  Consequently, 
a  certain  amount  of  gas  bypasses  the  gas  flow  channels  through  the 
porous  media  (a.k.a.  cross-leakage  flow).  This  effect  may  lead  to  a 
higher  concentration  of  the  reactants  in  the  regions  close  to  the  RL 
and  might  be  a  reason  for  higher  current  density  values  at  these 
locations  (see  Figs.  7  and  8). 

For  a  25  A  load  current,  the  net  pressure  loss  was  700  Pa  over  the 
cathode  gas  flow  channel  and  80  Pa  over  the  anode  gas  flow  chan¬ 
nel.  These  pressure  losses  are  low  compared  to  similar  LTPEMFC 
operating  conditions;  nevertheless,  it  must  be  noted  that  although 
the  reported  values  are  the  net  pressure  losses  over  the  exact  length 
of  the  flow-field.  The  simulated  values  are  lower  than  the  measured 
ones  because  this  FITPEMFC  model  does  not  account  for  any  addi¬ 
tional  peripheral  pressure  losses.  For  the  same  load  current  (25  A), 
the  measured  pressure  loss  was  2400  Pa  at  the  cathode  side  and 
250  Pa  at  the  anode  side.  In  fact,  the  presence  of  the  gas  connectors, 
the  gas  pipes  towards  the  differential  pressure  transmitters,  the  in- 
,  and  outlet  manifold,  as  well  as  the  gas  distributor  lead  to  higher 
pressure  losses  as  gas  recirculation  zones,  expansion  zones  (e.g. 
sudden  expansion),  and/or  contraction  zones  may  exist.  Results 
elucidate  that  a  detailed  modeling  of  the  entire  gas  piping  system 


would  give  a  clearer  picture  of  absolute  pressure  values.  Over¬ 
all,  qualitative  results  were  in  agreement  with  the  particle  image 
velocimetry  (PIV)  measurements  presented  in  [32]. 

5.3.  Current  density  distribution 

Fig.  7  shows  the  current  density  distribution  at  the  MEA  surface 
(RL)  for  a  cell  operated  at  base  case  operating  conditions  at  5  A,  1 0  A, 
15  A,  20  A,  and  25  A  load  currents  (left  to  right).  All  five  operating 
conditions  have  similar  current  density  distributions  with  higher 
values  towards  the  air  inlet  and  lower  values  towards  the  air  out¬ 
let.  Oxygen  decreases  almost  linearly  along  the  gas  flow  channel 
and  towards  the  RL  due  to  electrochemical  reactions.  The  distri¬ 
bution  of  the  oxygen  mass  fraction  dominates  the  current  density 
distribution.  Additionally,  the  distribution  is  strongly  influenced 
by  the  bipolar  plate  or  flow-field  structure.  At  the  exit  of  the  cath¬ 
ode  gas  flow  channel,  the  O2  mole  fraction  is  calculated  to  be  0.15 
for  the  given  cathode  stoichiometries  and  inlet  mole  fractions.  It 
is  clear  that  no  water  flooding  occurs  within  the  porous  media, 
and  the  water  vapor  mass  fraction  increases  towards  the  cathode- 
side  outlet  (oxygen  dilution).  Fligher  values  are  observed  under  the 
land  areas  and  close  to  the  borders  of  the  MEA,  especially  towards 
the  cathode  outlet  region.  This  is  reasonable  because  the  water 
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Fig.  7.  Cathode-side  RL  current  density  distribution  in  the  x-y-plane  (z  =  -379  x  10-6  m)  at  base  case  operating  conditions  for  5  A,  10  A,  15  A,  20  A,  and  25  A  load  currents 
(left  to  right). 
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Table  5 

Mean  solid-phase  temperature  within  the  HTPEMFC  components  /  °C. 


1  A 

5  A 

10  A 

15  A 

20  A 

25A 

Anode  side 

Cu 

160.07 

160.19 

160.329 

160.442 

160.556 

160.64 

BPP 

160.04 

160.193 

160.35 

160.491 

160.614 

160.73 

GDL 

160.02 

160.5 

161.02 

161.47 

161.91 

162.25 

RL 

160.06 

160.69 

161.39 

161.99 

162.58 

163.03 

MEM 

160.079 

160.75 

161.48 

162.126 

162.75 

163.22 

Cathode  side 
RL 

160.08 

160.791 

161.55 

162.226 

162.87 

163.41 

GDL 

160.055 

160.57 

161.13 

161.63 

162.11 

162.47 

BPP 

160.04 

160.626 

160.78 

160.924 

161.048 

161.16 

Cu 

160.07 

160.623 

160.76 

160.87 

160.99 

161.07 

vapor  has  to  diffuse  through  the  porous  media  towards  the  gas 
flow  channel  before  it  leaves  the  cell. 

In  addition  to  the  oxygen  mass  fraction,  the  current  density  is 
influenced  by  the  fluid-phase  (and/or  solid-phase)  temperature  as 
cold  gases  enter  the  cell  (compare  with  Figs.  9  and  10).  It  can  be 
seen  that  the  influence  is  more  significant  at  higher  currents  due  to 
larger  total  gas-flow  rates.  This  effect  is  even  noticed  in  the  region 
of  the  anode  inlet  for  higher  load  currents.  Over  the  entire  MEA 
area,  higher  current  density  values  are  observed  under  the  ribs  than 
under  the  channels  of  the  bipolar  plate.  The  maximum  current  den¬ 
sity  values  are  seen  in  the  regions  where  the  1 80°  bends  are  located. 
This  may  be  caused  by  higher  local  velocity  and  reactant  concen¬ 
tration  values  within  the  RL,  as  discussed  above  and  mentioned  in 
[23,32]. 

Fig.  8  provides  fine  details  about  the  current  density  distribution 
along  the  x-axis  and  highlights  the  above-mentioned  aspects.  Close 
to  the  cathode  inlet,  the  current  density  under  the  land  is  much 
higher  than  under  the  channel.  Moreover,  the  observed  variations 
become  less  pronounced  towards  the  cathode  outlet  for  all  load 
currents.  The  individual  current  density  peaks  at  each  180°  bend 
can  clearly  be  seen.  The  above  results  indicate  that  the  channel- 
to-land  ratio  should  be  seen  as  an  optimization  parameter  for  the 
flow-field  layout  used  for  FITPEMFC,  especially  at  the  cathode  side. 

As  for  the  anode  side,  it  was  observed  that  the  distribution  of 
the  hydrogen  mass/mole  fraction  is  much  smoother,  having  higher 
values  in  the  gas  flow  channels  than  in  the  porous  media  (con¬ 
sumption  in  z-direction  towards  the  RL).  The  influence  of  the  BPP 
structure  can  hardly  be  seen.  It  is  expected  that  hydrogen  easily 
reaches  the  catalysts,  where  it  is  finally  consumed.  At  the  exit  of 
the  anode  gas-flow  channel,  the  H2  mole  fraction  is  calculated  to 
be  0.985  for  the  given  anode  stoichiometries  and  inlet  mole  frac¬ 
tions.  This  variation  is  small  due  to  the  low  molar  mass  of  H2  and, 
because  no  water  is  produced  at  the  anode  side.  Additionally,  it  is 
assumed  that  no  water  is  transferred  from  the  cathode  to  the  anode 
side  via  the  membrane  or  water  absorbed  by  the  membrane. 

5.4.  Solid-phase  temperature  distribution 

The  solid-phase  temperature  within  the  gold-plated  copper  cur¬ 
rent  collectors  and  both  BPPs  seems  to  be  almost  uniform  due  to  the 
high  thermal  conductivities  of  both  materials  compared  to  all  other 
thermal  conductivities  within  the  setup.  Only  slightly  higher  solid- 
phase  temperatures  are  observed  within  these  components  when 
drawing  higher  load  currents.  Table  5  lists  the  calculated  values 
for  the  anode-  and  cathode-side  components.  The  mean  RL  solid- 
phase  temperature  rise  at  25  A  load  current  at  the  cathode  side  is 
3.41  °C  and  2.87  °C  at  20  A,  compared  to  no  load  operating  condi¬ 
tions  (set  to  160  °C).  It  is  expected  that  this  rise  in  temperature  will 
be  much  more  pronounced  at  load  currents  of  30-40  A,  since  much 
more  head  will  produced.  Anyways,  these  operating  points  are  not 
explicitly  discussed  in  this  work. 


Within  the  porous  media,  the  solid-phase  temperature  is 
strongly  influenced  by  the  cold  gases  entering  the  cell  (Fig.  9).  This 
fact  is  observed  at  both  gas  inlets.  At  the  anode  side,  even  though  the 
total  gas  flow  is  lower,  the  effect  is  much  more  pronounced  in  terms 
of  temperature  than  on  the  cathode  side.  On  the  other  hand,  a  much 
larger  MEA  area  is  influenced  on  the  cathode  side.  The  maximum 
solid-phase  temperature  is  observed  close  to  the  region  where  the 
highest  current  density  is  also  located  (164  °C  at  a  20  A  load  cur¬ 
rent  and  1 65  °C  at  a  25  A  load  current).  When  taking  a  closer  look  at 
the  solid-phase  temperature  distribution,  it  is  seen  that  the  shape 
of  both  flow-fields  influences  its  distribution.  The  temperature  in 
the  regions  under  the  gas-flow  channels  is  somewhat  lower.  This 
temperature  difference  almost  vanishes  towards  the  exit.  For  base 
case  operating  conditions  and  a  20  A  load  current,  the  total  amount 
of  heat  produced  is  calculated  to  be  11.22  W.  Most  of  the  heat  is 
produced  within  the  cathode-side  RL  (10.61  W).  Of  this  value,  the 
irreversible  reaction  heat  comprises  69%,  the  reaction  entropy  com¬ 
prises  30%,  and  the  total  joule  heating  (resistive  heating)  makes  up 
the  remaining  1%. 

A  major  drawback  of  this  model  setup  is  the  fact  that  a  constant 
temperature  is  defined  at  both  copper  current  collector  bound¬ 
aries.  In  fact,  much  higher  temperatures  are  expected  within  the 
endplates  (up  to  1 80-200  °C  in  the  region  close  to  the  heating  ele¬ 
ments)  because,  in  most  single  cell  setups,  heating  elements  are 
used  to  keep  the  FITPEMFC  at  operating  temperature. 

5.5.  Fluid- (gas )-phase  temperature  distribution 

The  gases  are  heated  as  they  diffuse  through  the  porous  media 
(in  the  z-direction  towards  the  RL).  The  fluid-(gas)-phase  tempera¬ 
ture  distribution  for  different  load  currents  within  the  cathode-side 
RL  is  given  in  Fig.  10.  Again,  when  drawing  more  current,  more  gas 
enters  the  cell.  Consequently,  it  takes  longer  until  the  gas  reaches 
the  cell  operating  temperature.  More  heat  needs  to  be  exchanged 
between  the  two  phases,  and  the  thermal  equilibrium  between  the 
solid-,  and  the  fluid-(gas)-phases  shifts  in  the  direction  of  the  chan¬ 
nel.  Within  the  anode-side  porous  media,  0.448  W  are  transferred, 
whereas  1.619W  are  transferred  between  the  two  phases  at  the 
cathode-side. 

Fig.  11  shows  the  anode  and  cathode-side  fluid-(gas)-phase 
temperature  behavior  within  the  middle  of  the  gas  flow  channel. 
At  the  anode  side,  hydrogen  heats  up  very  quickly  (high  thermal 
conductivity  and  a  high  heat  capacity)  and  reaches  cell  operating 
temperature  shortly  after  entering  the  gas  flow  channels.  At  the 
cathode  side,  it  takes  much  longer  for  the  gas  mixture  to  reach 
thermal  phase  equilibrium.  Both,  Fig.  11  (z  =  500  x  10-6  m)  and 
Fig.  10  (z=  -379  x  10-6  m)  show  the  same  trend  when  it  comes  to 
fluid-phase  temperature  distribution.  The  difference  is  that  it  takes 
slightly  longer  until  phase  equilibrium  is  reaches  within  the  gas 
flow  channels  (higher  convection)  than  within  the  porous  media. 
This  behavior  was  observed  for  all  load  currents. 

The  heat  flux  vectors  within  the  high  magnification  image  in 
Fig.  1 1  depict  the  solid-to-fluid  heat  transfer  via  the  gas  flow  chan¬ 
nel  (bipolar  plate)  walls  (using  a  constant  heat  transfer  coefficient 
at  these  boundaries).  It  can  be  seen  that  the  gas  temperature  is 
becomes  slightly  higher  close  to  the  bipolar  plate  (gas  flow  channel) 
boundaries.  For  base  case  operating  conditions  and  a  20  A  load  cur¬ 
rent,  0.01 46  W  and  1.0747W  are  transferred  through  these  walls 
at  the  anode  and  cathode  side,  respectively. 

5.6.  Membrane  conductivity 

Eqs.  ( 1 1 )-( 1 3 )  satisfactorily  describe  the  overall  membrane  con¬ 
ductivity  only  in  a  small  temperature  range  around  160°C.  For 
higher  temperatures,  e.g.,  TS>170°C,  the  calculated  conductivi¬ 
ties  do  not  match  the  measured  data  in  [34],  indicating  that,  in 
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Fig.  9.  Cathode-side  RL  solid-phase  temperature  distribution  in  x-y-plane  (z  =  -379  x  1 0-6  m)  at  base  case  operating  conditions  for  5  A,  1 0  A,  1 5  A,  20  A,  and  25  A  load  current 
(left  to  right). 


Fig.  10.  Cathode-side  RL  fluid-(gas)-phase  temperature  distribution  inx-y-plane  (z=-379  x  10-6  m)  at  base  case  operating  conditions  for  5  A,  10  A,  15  A,  20,  and  25  A  load 
currents  (left  to  right). 


this  range  of  temperatures,  other  mechanisms  of  proton  conduc¬ 
tion  may  be  present  and  may  be  influenced  by  the  concentration 
changes  of  the  acid  within  the  PBI/H3PO4  system.  Though  the  con¬ 
ductivity  values  of  the  PBI/H3PO4  sol-gel  membranes  are  quite 


high,  the  values  are  still  lower  than  the  reported  values  for  hot 
concentrated  H3PO4  by  a  factor  of  approximately  three  (160  °C). 
This  indicates  that  the  transport  process  through  this  type  of  mem¬ 
brane  is  not  as  continuous  as  in  hot  concentrated  H3P04,  which  is 
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Fig.  11.  Fluid-(gas)-phase  temperature  distribution  within  the  anode-  and  cathode-side  gas  flow  channels  (z  =  500  x  1 0-6  m)  at  base  case  operating  conditions  in  the  x-y-plane 
for  5  A  (anode-side),  20A  (anode-side),  5  A  (cathode-side),  and  20  A  (cathode-side)  load  currents  (left  to  right). 


C.  Siegel  et  al.  /  Journal  of  Power  Sources  196(2011)  2735-2749 


2747 


a  in 


c  out 


a  in  c  out 

2- 10-3 
1.98- 10-3 
1.96-10-3 
1.94-10-3 
1.92-10-3 
1.90-10-3 
1.88-10-3 
1.86-10-3 
1.84-10-3 
1.82-10-3 
1.80-10-3 


•  -  16.4 
-  16.2 
-  16 

15.8 

•  -  15.6 
-  15.4 
I  15.2 
•  15 


c  in  a  out  c  in  a  out 

Fig.  12.  (a)  Membrane  conductivity  at  various  solid-phase  temperatures;  (b)  membrane  resistance  distribution  in  the  x-y-plane  (z=  -430  x  10-6  m)  for  base  case  operating 
conditions  at  a  20  A  load  current;  (c)  membrane  conductivity  distribution  in  the  x-y-plane  (z  =  -430  x  10-6  m)  for  base  case  operating  conditions  at  a  20  A  load  current. 


consistent  with  other  types  of  high-temperature  stable  membranes 
[7]. 

This  model  assumes  no  gas  crossover  or  water  transfer  through 
the  membrane.  In  Eq.  (12),  the  pre-exponential  factor  and  the  acti¬ 
vation  energy  are  taken  to  be  constant.  Consequently,  the  local 
membrane  conductivity  distribution  follows  the  solid-phase  tem¬ 
perature  distribution.  The  minimum  (15.7Sm-1)  is  close  to  the 
anode-side  inlet,  whereas  the  maximum  (17  Sm-1)  is  located  in 
the  lower  third  of  the  MEA,  as  can  be  seen  from  Fig.  12. 

The  membrane  resistance  calculated  from  the  model  was 
1.877  m£2,  which  is  in  good  agreement  with  the  measured  value. 
Again,  its  distribution  follows  the  membrane  conductivity  distri¬ 
bution.  At  base  case  operating  conditions  and  a  20  A  load  current, 
the  voltage  loss  over  the  membrane  was  0.03927  V. 


5.7.  Changes  in  H3PO4  during  cell  operation 

This  model  assumes  no  water  through  or  water  uptake  by  the 
membrane  (see  assumption  6).  It  uses  a  simple  approach  to  couple 
the  partial  pressure  of  water  to  the  phosphoric  acid  concentration 
and  temperature  within  both  RL,  as  seen  in  Eq.  (25).  The  calcu¬ 
lated  H3PO4  (and  amorphous  phase  H3PO4)  concentration  changes 
for  different  load  currents  because  more  water  is  produced.  At  the 
anode  side,  the  partial  pressure  of  water  remains  more  or  less  con¬ 
stant,  whereas  it  increases  towards  the  exit  at  the  cathode  side.  The 
model  results  indicate  that  the  H3PO4  concentration  is  reduced  for 
higher  water  vapor  partial  pressures.  For  base  case  operating  con¬ 
ditions,  it  is  observed  that  the  H3P04  concentration  reaches  much 
higher  values  than  the  initially  assumed  concentration  of  85  wt.% 
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Fig.  13.  Phosphoric  acid  concentration  distribution  in  the  x-y-plane  (z=  -379  x  10-6  m)  for  base  case  operating  conditions  at  5  A  (left)  and  20  A  load  currents  (right). 


at  RTP,  as  can  be  seen  in  Fig.  13.  At  the  anode  and  cathode  side, 
the  mean  concentrations  are  calculated  to  be  100.23%  and  96.79%, 
respectively  (based  on  molecular  H3P04).  The  mean  partial  pres¬ 
sure  of  water  is  calculated  to  be  1802  Pa  at  the  anode-side  and 
18134  Pa  at  the  cathode-side  RL  (base  case  operating  conditions 
and  a  20  A  load  current).  Moreover,  from  Fig.  13  it  can  be  seen 
that  the  highest  values  at  the  cathode  side  are  located  close  to 
the  inlet  in  the  lower  third  over  the  MEA.  Locally,  lower  values  are 
seen  under  the  land  areas  because  more  water  vapor  is  present.  It 
must  be  noted  that  the  reported  concentration  values  may  change 
when  accounting  for  water  adsorption  and  desorption  in  PBI/H3PO4 
systems. 

The  FI3PO4  concentration  values  and  the  temperature  influence 
the  gas  diffusivity  and  gas  solubility  in  amorphous  phase  H3  PO4.  For 
base  case  operating  conditions,  slightly  lower  air  diffusivity  values 
were  located  close  to  the  cathode  inlet,  whereas  almost  constant 
values  were  observed  over  the  remaining  MEA  area.  All  calculated 
diffusivity  values  were  in  good  agreement  with  the  values  reported 
for  hot  concentrated  H3PO4. 

The  air  solubility  had  the  same  characteristic  distribution  as  the 
phosphoric  acid  concentration  in  Fig.  13.  For  base  case  operating 
conditions  and  a  20  A  load  current,  the  solubility  ranged  from  0.504 
to  0.529  mol  m-3  atm-1.  This  means  that  both  the  gas  diffusivity 
and  the  solubility  are  somewhat  lower  close  to  the  cathode  inlet, 
influencing  the  local  current  density  values. 

Nevertheless,  much  more  experimental  and  theoretical  work 
(e.g.,  molecular  dynamics  simulations)  is  necessary  to  precisely 
identify  the  interactions  of  FI3PO4,  PBI,  and  water  vapor  to  develop 
an  adequate  phosphoric  acid  transport  model  for  continuous  oper¬ 
ation  and  start/stop  cycling. 

6.  Conclusion 

A  complete  non-isothermal  three-dimensional  model  of  a 
FITPEMFC  setup  using  a  PBI/FI3PO4  sol-gel  membrane  was 


developed,  modeled,  and  solved  using  COMSOL  Multiphysics®. 
Additional  equations  were  directly  coded  into  the  FEM/CFD 
software  using  scalar,  boundary,  and  subdomain  expressions.  Com¬ 
putational  aspects  are  listed  in  order  to  paint  a  clear  picture  of 
the  soft-,  and  hardware  requirements  when  solving  large  scale 
models  using  finite  elements.  Electrochemical  reactions  were 
described  using  an  agglomerate  approach  that  considered  the 
diffusivity  and  solubility  of  the  gases.  The  model  was  able  to 
reproduce  the  achieved  experimental  results  for  0-20  A  load  cur¬ 
rents.  The  membrane  conductivity  was  modeled  using  an  Arrhenius 
approach  that  seems  to  be  valid  for  solid-phase  temperatures  up 
to  150-1 60  °C.  For  higher  temperatures,  this  approach  may  over¬ 
predict  the  membrane  conductivity  since  the  amorphous  phase 
FI3PO4  concentration  changes  may  influence  the  proton  conduc¬ 
tion  mechanisms.  Results  also  show  how  the  FI3PO4  concentration 
is  influenced  during  cell  operation.  Beside,  the  model  setup  is  able 
to  simulate  the  dependency  of  the  current  density  distribution  on 
the  fluid-flow  distribution  and  on  the  solid-  and  fluid-(gas)-phase 
temperature  distribution  in  a  FITPEMFC.  Flighest  solid-phase  tem¬ 
perature  is  observed  in  the  region  of  the  highest  current  density. 
Moreover,  it  highlights  the  interaction  between  the  two  tempera¬ 
tures,  including  localized  cooling  effects  close  to  the  gas  inlets  by 
reporting  values  for  the  amount  of  transferred  heat.  Consequently, 
the  fluid-(gas)-phase  inlet  temperature  and  the  gas  manifold  and 
flow-field  design  should  be  seen  as  an  important  factor  for  optimal 
operation,  especially  when  designing  gas  manifolds  and  gas-inlet 
sections  for  larger  FITPEMFC  stacks.  This  general-purpose  model 
may  be  useful  when  analyzing  FITPEMFC  operational  behavior 
because  it  exemplifies  the  complex  couplings  between  all  trans¬ 
port  equations  and  electrochemical  reactions.  Based  on  the  results, 
it  is  now  possible  to  optimize  the  flow-field  structure  and  minimize 
current  density  distribution  gradients  over  the  MEA.  It  will  be  the 
basis  of  future  FITPEMFC  modeling  activities  that  may  include  pos¬ 
sible  water  transfer  through  the  PBI/FI3PO4  system,  gas  crossover 
and  an  ameliorated  FI3PO4  transport  model. 
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